#Use renv to install the requisite packages for this replication.
#install.packages("renv")
#renv::install()
source("preamble.R")

IMG_FILE_EXT = ".pdf"

# Completion Rate by Treatment Condition ----------------------------------
message("Completion Rate by Treatment Condition")

t1 <- combined %>%
  tabyl(treatment, status) %>%
  flextable() %>%
  set_header_labels(
    treatment = "Treatment\nGroup",
    complete = "Number of\nCompletes",
    dropout = "Number of\nDropouts"
  ) %>%
  theme_vanilla() %>% 
  bg(
    bg = "white",
    part = "all"
  ) %>% 
  set_table_properties(layout = "autofit", width = 1) %>%
  set_caption(
    number_tab("Unit Completion Status by Treatment")
  )

save_as_image(x = t1, "figures/a_tab_complete_by_treat.png")

# Interview Duration by Treatment Condition --------------------------------
message("Interview Duration by Treatment Condition")

t2 <- combined %>%
  group_by(treatment) %>%
  mutate(duration = duration/60) %>%
  summarize(
    mean = mean(duration, na.rm = TRUE),
    median = median(duration, na.rm = TRUE),
    min = min(duration, na.rm = TRUE),
    max = max(duration, na.rm = TRUE),
    q01 = quantile(duration, 0.01, na.rm = TRUE),
    q25 = quantile(duration, 0.25, na.rm = TRUE),
    q75 = quantile(duration, 0.75, na.rm = TRUE),
    q99 = quantile(duration, 0.99, na.rm = TRUE)
  ) %>%
  adorn_rounding(digits = 2) %>% 
  flextable() %>% 
  set_header_labels(
    treatment = "Treatment\nGroup",
    mean = "Mean",
    median = "Median",
    min = "Minimum",
    max = "Maximum",
    q01 = "1st Percentile",
    q25 = "25th Percentile",
    q75 = "75th Percentile",
    q99 = "99th Percentile"
  ) %>% 
  theme_vanilla() %>% 
  bg(
    bg = "white",
    part = "all"
  ) %>% 
  set_table_properties(
    layout = "autofit",
    width = 1
  ) %>% 
  set_caption(
    number_tab("Interview Duration by Treatment Condition (in minutes)")
  )

save_as_image(x=t2, "figures/a_tab_intv_dur_by_treat.png")

# Completion Time by Treatment Condition ----------------------------------
message("Completion Time by Treatment Condition")

t3 <- combined %>%
  filter(completed) %>%
  group_by(treatment) %>%
  mutate(duration = duration/60) %>%
  summarize(
    mean = mean(duration, na.rm = TRUE),
    median = median(duration, na.rm = TRUE),
    min = min(duration, na.rm = TRUE),
    max = max(duration, na.rm = TRUE),
    q01 = quantile(duration, 0.01, na.rm = TRUE),
    q25 = quantile(duration, 0.25, na.rm = TRUE),
    q75 = quantile(duration, 0.75, na.rm = TRUE),
    q99 = quantile(duration, 0.99, na.rm = TRUE)
  ) %>%
  adorn_rounding(digits = 2) %>% 
  flextable() %>% 
  set_header_labels(
    treatment = "Treatment\nGroup",
    mean = "Mean",
    median = "Median",
    min = "Minimum",
    max = "Maximum",
    q01 = "1st Percentile",
    q25 = "25th Percentile",
    q75 = "75th Percentile",
    q99 = "99th Percentile"
  ) %>% 
  theme_vanilla() %>% 
  bg(
    bg = "white",
    part = "all"
  ) %>% 
  set_table_properties(
    layout = "autofit",
    width = 1
  ) %>% 
  set_caption(
    number_tab("Completion Time by Treatment Condition (in minutes)")
  )

save_as_image(x=t3, "figures/a_tab_compl_time_by_treat.png")

# Composition -------------------------------------------------
message("Composition")

add_pcts <- function(df) {
  df %>% group_by(treatment) %>% mutate(pct = n/sum(n)) %>% ungroup()
}

d <- bind_rows(combined %>%
                 mutate(treatment = case_when(
                   treatment == "Conf." ~ "Treatment 2\n(Conf.)",
                   treatment == "Elab/Rel." ~ "Treatment 1\n(Elab/Rel.)",
                   TRUE ~ "Control"
                 )),
               combined %>% 
                 mutate(treatment = "Overall"))

bind_rows(
  d %>%
    count(treatment, val=treatment, var="Condition") %>%
    add_pcts() %>%
    arrange(val != "Overall", val != "Control") %>%
    mutate(val = gsub("\n.*", "", val)),
  d %>%
    mutate(device_type = forcats::fct_relevel(device_type, 
                                              "desktop", "smartphone","tablet", "other")) %>%
    count(treatment, val=device_type, var="Device") %>%
    add_pcts(),
  d %>%
    count(treatment, val=gender, var="Gender") %>%
    add_pcts(),
  d %>%
    count(treatment, val=age5, var="Age") %>%
    add_pcts(),
  d %>%
    count(treatment, val=educ5, var="Educ") %>%
    add_pcts(),
  d %>%
    count(treatment, val=hh_inc5, var="Household\nIncome") %>%
    add_pcts(),
  d %>%
    count(treatment, val=work_status, var="Employed") %>%
    add_pcts()
) %>%
  mutate(val = val %>%
           as.character() %>%
           replace_na("Refused") %>%
           stringr::str_to_title() %>%
           stringr::str_replace("Hs", "HS")) %>%
  arrange(val == "Refused", val == "No") %>%
  mutate(val = as_factor(val), 
         var = as_factor(var),
         treatment = factor(treatment, 
                            levels = c("Overall", 
                                       "Control", 
                                       "Treatment 1\n(Elab/Rel.)",
                                       "Treatment 2\n(Conf.)"))) %>%
  ggplot(aes(y=val, x=treatment, fill=pct*100)) +
  geom_tile(color="white") +
  geom_text(aes(label = sprintf("n=%d", n),
                color = ifelse(pct > 0.5, "black", "white")),
            fontface = "bold", size = 4) +
  scale_x_discrete(name = "Condition") +
  scale_y_discrete(limits = rev, name = "") +
  scale_color_identity() +
  scale_fill_viridis_c(name = "% in Group:") +
  facet_grid(var ~ ., scales= "free_y", space = "free_y", labeller=label_wrap_gen(width = 10)) +
  theme_bw() +
  theme(legend.title = element_text(size=12, vjust=0.8),
        legend.position = "bottom",
        strip.text.x = element_text(size=22, lineheight=1),
        panel.grid = element_blank(),
        strip.text.y = element_text(size=17, hjust=0, lineheight=1, angle=0),
        axis.text = element_text(size=16),
        axis.title = element_text(size=18, lineheight=1))

ggsave(here(paste0("figures/a_composition", IMG_FILE_EXT)), 
       height=12, width=10)


# Summary of Human-Coded Quality Criteria ---------------------------------
message("Summary of Human-Coded Quality Criteria")

coded_qual_resolved %>%
  group_by(q, treat) %>%
  summarise(across(all_of(quality_criteria), ~mean(.x, na.rm=T)), 
            .groups = "drop") %>%
  pivot_longer(cols=all_of(quality_criteria)) %>%
  mutate(value = ifelse(q %in% c("news", "occu") & name == "explanatory", NA, value)) %>%
  mutate(treat = fct_recode(treat, "Elab/Rel."="Elaboration", "Conf."="Confirmation"),
         name = stringr::str_to_title(name),
         q = exp_ordered_labels_f(q)) %>%
  mutate(name = factor(name, levels=stringr::str_to_title(quality_criteria))) %>%
  ggplot(aes(x=treat, y=value, fill=treat)) +
  facet_grid(name ~ q, scales = "free_x", labeller=label_wrap_gen(width = 10)) +
  labs(caption = str_wrap(paste("Note: News",
                                "and occupation experiments did not have",
                                "an open-ended control question for comparison."), 120),
       x = "Treatment Condition") +
  geom_bar(stat="identity") +
  geom_text(aes(label=ifelse(is.na(value), "N/A", sprintf("%0.1f%%", 100*value)),
                y=ifelse(is.na(value), 0.25, value), vjust=ifelse(value > 0.5, 1.5, -.25)), 
            color="black", size=4) +
  coord_cartesian(clip = "off") +
  scale_y_continuous(labels = scales::percent_format(1),
                     name = "Rate of Incidence") +
  scale_fill_manual(values=c("grey","grey")) +
  theme_bw() +
  theme(legend.position = "none",
        panel.background = element_rect(),
        plot.caption = element_text(lineheight = 1, size = 10, hjust = 0),
        strip.text = element_text(size=22, lineheight=1),
        strip.text.y = element_text(angle=0, hjust = 0),
        axis.text.y = element_text(size=14),
        axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
        axis.title = element_text(size=18, lineheight=1),
        panel.spacing = unit(1, "lines"),
        plot.margin = margin(t = 10, r = 10, b = 10, l = 10))

ggsave(here(paste0("figures/a_summary_qual_criteria", IMG_FILE_EXT)), 
       height=12, width=12)

# Summary of NLP-Based Informational Measures -----------------------------
message("Summary of NLP-Based Informational Measures")

nlp_measures_l <- nlp_measures %>%
  mutate(response = case_when(
    treatment == "Control" & grepl("econd_reason_1", response) ~ "econd_last_response",
    TRUE ~ response
  )) %>%
  mutate(exp = gsub("_(combined|last)_response", "", response) %>%
           exp_ordered_labels_f()) %>%
  pivot_longer(names_to = "measure",
               values_to = "val",
               unique_words:shannon_entropy)

nlp_measures_l %>%
  filter(grepl("last_response", response) & !(grepl("why_frustrated", response))) %>%
  filter(!(grepl("(News|Occu)", exp) & grepl("Control", treatment))) %>%
  mutate(measure = nlp_measures_ordered_labels_f(measure)) %>%
  ggplot(aes(x=treatment, y=val)) +
  geom_boxplot() +
  labs(caption=str_wrap(paste("Note: Measures",
                        "are shown for the post-probing response text in the elab/rel. treatment condition"), 150),
       x="Treatment Condition"
       ) +
  scale_y_continuous(name = "Distribution of Measures") +
  scale_fill_manual(values=c("grey", "grey")) +
  facet_grid(measure ~ exp, 
             scales = "free",
             , labeller=label_wrap_gen(width = 10)) +
  theme_bw() +
  theme(strip.text = element_text(size=17),
        legend.position = "bottom",
        plot.caption = element_text(lineheight = 1, size = 10, hjust = 0),
        panel.background = element_rect(),
        strip.text.y = element_text(size=18, angle=0, hjust=0, lineheight=1),
        axis.text.x = element_text(size=16, angle=40, hjust=1),
        axis.text = element_text(size=18),
        axis.title = element_text(size=16, lineheight=1))

ggsave(here(paste0("figures/a_summary_nlp_measures", IMG_FILE_EXT)), 
       height=12, width=12)

# Summary of Respondent Attrition -----------------------------------------
message("Summary of Respondent Attrition")

combined %>%
  select(treatment, matches("rightafter")) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  pivot_longer(-treatment) %>%
  mutate(exp = name %>%
           str_remove("_dropout.*") %>%
           exp_ordered_labels_f() %>% 
           factor(., levels = unique(.))) %>%
  ggplot(aes(x=treatment, y=value)) +
  facet_grid(exp ~ ., labeller=label_wrap_gen(width = 10)) +
  scale_x_discrete(name = "Treatment Group") +
  scale_y_continuous(name = "% of Sample Attrition Out Right After Question", 
                     labels = scales::percent_format(1)) +
  geom_bar(stat = "identity", fill = "grey") +
  geom_text(aes(label=sprintf("%0.1f%%", 100*value),
                y=value, vjust=ifelse(value > 0.04, 1.5, -.25)), 
            color="black", size=5) +
  theme_bw() +
  theme(panel.background = element_rect(),
        strip.text = element_text(size=22, lineheight = 1),
        strip.text.y = element_text(size=18, angle=0, hjust=0.1),
        axis.text = element_text(size=14),
        axis.title = element_text(size=18, lineheight=1),
        panel.spacing = unit(1, "lines"),
        plot.margin = margin(t = 10, r = 10, b = 10, l = 10))

ggsave(here(paste0("figures/a_summary_attrition", IMG_FILE_EXT)), 
       height=8, width=10)

# Summary of Self-Reported User Experience --------------------------------
message("Summary of Self-Reported User Experience")

add_pcts <- function(df) {
  df %>% 
    group_by(treatment) %>% 
    mutate(pct = n/sum(n)) %>% 
    ungroup()
}

d <- bind_rows(
  combined,
  combined %>% mutate(treatment = "Overall")
)

bind_rows(
  d %>%
    count(treatment, val=quality, var="Quality", .drop = FALSE) %>%
    mutate(var = "Quality") %>%
    add_pcts(),
  # d %>%
  #   count(treatment, val=formal_tone, var="Formal Tone", .drop = FALSE) %>%
  #   add_pcts(),
  d %>%
    count(treatment, val=ease, var="Ease", .drop = FALSE) %>%
    mutate(var = "Ease") %>%
    add_pcts(),
  d %>%
    count(treatment, val=frustration, var="Frustration", .drop = FALSE) %>%
    mutate(var = "Frustration") %>%
    add_pcts(),
  d %>%
    count(treatment, val=satisfaction, var="Satisfaction", .drop = FALSE) %>%
    mutate(var = "Satisfaction") %>%
    add_pcts()
) %>%
  mutate(val = val %>%
           as.character() %>%
           replace_na("Missing") %>%
           stringr::str_to_title() %>%
           stringr::str_replace("Hs", "HS")) %>%
  arrange(val == "Missing") %>%
  mutate(val = as_factor(val), 
         var = as_factor(var),
         treatment = factor(treatment, 
                            levels = c("Overall", "Control", "Conf.", "Elab/Rel."))) %>%
  ggplot(aes(y=val, x=treatment, fill=pct*100)) +
  geom_tile(color="white") +
  geom_text(aes(label = sprintf("n=%d", n),
                color = ifelse(pct > 0.5, "black", "white")),
            fontface = "bold", size = 5) +
  scale_x_discrete(name = "Group") +
  scale_y_discrete(limits = rev, name = "") +
  scale_color_identity() +
  scale_fill_viridis_c(name = "% in Group:") +
  facet_grid(var ~ ., scales= "free_y", space = "free_y") +
  theme_bw() +
  theme(legend.title = element_text(size=12, vjust=0.8),
        legend.text = element_text(size=12),
        legend.position = "bottom",
        strip.text.x = element_text(size=22, lineheight=1),
        panel.grid = element_blank(),
        strip.text.y = element_text(size=17, lineheight=1),
        axis.text = element_text(size=16),
        axis.title = element_text(size=18, lineheight=1))

ggsave(here(paste0("figures/a_summary_user_experience", IMG_FILE_EXT)), 
       height=12, width=12)

# Detection Accuracy by Device Type ---------------------------------------
message("Detection Accuracy by Device Type")

acc_dt <- map_dfr(unique(treat_conf$device_type), function(dt) {
  treat_conf_dt <- filter(treat_conf, device_type == dt)
  data.frame(
    device.type = paste0(dt, " (n=", nrow(treat_conf_dt), ")"),
    issue = with(treat_conf_dt, mean(issue_detect_correct, na.rm = T)),
    # econd_tone.acc = with(treat_conf, mean(econd_tone_detect_correct, na.rm = T)),
    # excluding when we just guess negative tone
    econd_tone = with(treat_conf_dt, mean(econd_tone_detect_correct[econd_tone_1_targets_detected_n == 1], na.rm = T)),
    econd_reason = with(treat_conf_dt, mean(econd_reason_detect_correct, na.rm = T)),
    news = with(treat_conf_dt, mean(news_detect_correct, na.rm = T)),
    occu = with(treat_conf_dt, mean(occu_detect_correct, na.rm = T))
  )
})
acc_dt %>% 
  mutate(device.type = stringr::str_to_sentence(device.type)) %>%
  flextable() %>%
  set_header_labels(
    values = append(list(device.type="Device",
                         metric="Detection\nPerformance\nMetric"),
                    exp_ordered_labels)
  ) %>%
  set_caption(
    number_tab("Detection Performance by Device")
  )

# Comparison of Response Categories (Detected in Treatment 2 vs. Answered in Control Close-End) ----
message("Comparison of Response Categories (Detected in Treatment 2 vs. Answered in Control Close-End)")

ggarrange(
  make_category_correlation_plot(cat_distr$econd_tone, 
                                 "Detected in Conf.\nCondition OE", 
                                 "Answered by Respondent\nin Control Condition CE") + 
    ggrepel::geom_text_repel(aes(label=cat)) +
    ggtitle("Econ. Cond. (Sent.)"),
  make_category_correlation_plot(cat_distr$news, 
                                 "Detected in Conf.\nCondition OE", 
                                 "Answered by Respondent\nin Control Condition CE") + 
    ggrepel::geom_text_repel(aes(label=cat),
                             max.overlaps = 5) +
    ggtitle("Pref. News"),
  make_category_correlation_plot(cat_distr$occu, 
                                 "Detected in Conf.\nCondition OE",
                                 "Answered by Respondent\nin Control Condition CE") + 
    ggrepel::geom_text_repel(aes(label=cat),
                             box.padding = 1) +
    ggtitle("Main Occu.")
)
ggsave(here(paste0("figures/a_det_vs_closed", IMG_FILE_EXT)),
       height = 12, width = 12)


# Comparison of Response Categories (Detected in Treatment 2 vs. Detected in Other Conditions) ----
message("Comparison of Response Categories (Detected in Treatment 2 vs. Detected in Other Conditions)")

ggarrange(
  make_category_correlation_plot(cat_distr$issue, 
                                 "Detected in Conf.\nCondition OE", 
                                 "Detected in Control Condition OE") + 
    ggrepel::geom_text_repel(aes(label=cat),
                             max.overlaps = 5) +
    ggtitle("Most Imp. Issue"),
  make_category_correlation_plot(cat_distr$news, 
                                 "Detected in Conf.\nCondition OE", 
                                 "Detected in Elab/Rel. Condition OE") + 
    ggrepel::geom_text_repel(aes(label=cat),
                             max.overlaps = 5) +
    ggtitle("Pref. News"),
  make_category_correlation_plot(cat_distr$occu, 
                                 "Detected in Conf.\nCondition OE",
                                 "Detected in Elab/Rel. Condition OE") + 
    ggrepel::geom_text_repel(aes(label=cat),
                             box.padding = 1) +
    ggtitle("Main Occu.")
)
ggsave(here(paste0("figures/a_det_vs_det_other", IMG_FILE_EXT)),
       height = 12, width = 12)

# Comparison of Response Categories (Confirmed in Treatment 2 vs. Answered in Control Close-End) ----
message("Comparison of Response Categories (Confirmed in Treatment 2 vs. Answered in Control Close-End)")

ggarrange(
  make_category_correlation_plot(cat_distr$econd_tone, 
                                 "Confirmed by Respondent\nin Conf. Condition OE", 
                                 "Answered by Respondent\nin Control Condition CE") + 
    ggrepel::geom_text_repel(aes(label=cat)) +
    ggtitle("Econ. Cond. (Sent.)"),
  make_category_correlation_plot(cat_distr$news, 
                                 "Confirmed by Respondent\nin Conf. Condition OE", 
                                 "Answered by Respondent\nin Control Condition CE") + 
    ggrepel::geom_text_repel(aes(label=cat)) +
    ggtitle("Pref. News"),
  make_category_correlation_plot(cat_distr$occu, 
                                 "Confirmed by Respondent\nin Conf. Condition OE",
                                 "Answered by Respondent\nin Control Condition CE") + 
    ggrepel::geom_text_repel(aes(label=cat),
                             box.padding = 1) +
    ggtitle("Main Occu.")
)
ggsave(here(paste0("figures/a_conf_vs_close", IMG_FILE_EXT)),
       height = 12, width = 12)

# Comparison of Response Categories (Confirmed in Treatment 2 vs. Detected in Treatment 2) ----
message("Comparison of Response Categories (Confirmed in Treatment 2 vs. Detected in Treatment 2)")

ggarrange(
  make_category_correlation_plot(cat_distr$issue, 
                                 "Confirmed by Respondent\nin Conf. Condition OE", 
                                 "Detected in Conf.\nCondition OE") + 
    ggrepel::geom_text_repel(aes(label=cat)) +
    ggtitle("Most Imp. Issue"),
  make_category_correlation_plot(cat_distr$econd_tone, 
                                 "Confirmed by Respondent\nin Conf. Condition OE", 
                                 "Detected in Conf.\nCondition OE") + 
    ggrepel::geom_text_repel(aes(label=cat)) +
    ggtitle("Econ. Cond. (Sent.)"),
  make_category_correlation_plot(cat_distr$econd_reason, 
                                 "Confirmed by Respondent\nin Conf. Condition OE", 
                                 "Detected in Conf.\nCondition OE") + 
    ggrepel::geom_text_repel(aes(label=cat)) +
    ggtitle("Econ. Cond. (Reason)"),
  make_category_correlation_plot(cat_distr$news, 
                                 "Confirmed by Respondent\nin Conf. Condition OE", 
                                 "Detected in Conf.\nCondition OE") + 
    ggrepel::geom_text_repel(aes(label=cat), 
                             max.overlaps = 5) +
    ggtitle("Pref. News"),
  make_category_correlation_plot(cat_distr$occu, 
                                 "Confirmed by Respondent\nin Conf. Condition OE",
                                 "Detected in Conf.\nCondition OE") + 
    ggrepel::geom_text_repel(aes(label=cat),
                             box.padding = 1) +
    ggtitle("Main Occu.")
)
ggsave(here(paste0("figures/a_conf_vs_det", IMG_FILE_EXT)),
       height = 12, width = 16)

# Detection Confusion Matrices -----------------------------------------
message("Detection Confusion Matrices")

full_conf_matrices <- lapply(
  c("issue","econd_reason","news","occu"),
  function (exp) {
    treat_conf %>%
      rename_at(paste0(exp,"_detect_correct"), ~"correct") %>%
      rename_at(paste0(exp,"_confirm"), ~"conf") %>%
      rename_at(paste0(exp,"_to_confirm"), ~"asked") %>%
      rename_at(paste0(exp,"_confirm_other"), ~"other") %>%
      filter(!is.na(correct)) %>%
      mutate(correct = case_when(
        conf == "yes" ~ asked,
        TRUE ~ other
      )) %>%
      count(asked, correct, name = "count") %T>% {
        with(., stopifnot(all(!is.na(count))))
      } %>%
      group_by(asked) %>%
      mutate(pct = count/sum(count)) %>%
      ungroup()
  }) %>% 
  setNames(c("issue","econd_reason","news","occu"))

full_conf_matrices$issue %>%
  make_full_confusion_heatmap() +
  labs(#title="Confusion Matrix for Most Important Issue",
    caption="Note: Mutiple issues may have been detected for a response, in which case a single issue was sampled for probing",
    y="Issue Detected + Probed", x="Issue Confirmed")
ggsave("figures/a_heatmap_issue.png", height=10, width=12)

full_conf_matrices$econd_reason %>%
  make_full_confusion_heatmap() +
  labs(#title="Confusion Matrix for Economic Condition Reason",
    caption="Note: Mutiple reasons may have been detected for a response, in which case a single reason was sampled for probing",
    y="Reason Detected + Probed", x="Reason Confirmed") 
ggsave("figures/a_heatmap_econd.png", height=10, width=12)

full_conf_matrices$news %>%
  mutate(asked = fct_relevel(asked, "another_newspaper", "another_tv_network", after = Inf),
         correct = fct_relevel(correct, "another_newspaper", "another_tv_network", after = Inf)) %>%
  make_full_confusion_heatmap() +
  labs(#title="Confusion Matrix for Preferred News Source",
    caption="Note: Mutiple outlets may have been detected for a response, in which case a single outlet was sampled for probing",
    y="News Source Detected + Probed", x="News Source Confirmed") +
  theme(axis.text.x = element_text(angle=60))
ggsave("figures/a_heatmap_news.png", height=10, width=12)

full_conf_matrices$occu %>%
  # remove extra 's' at end of occupation labels  
  mutate( 
    asked = case_when(str_ends(asked, "s") ~ str_sub(asked, 1, -2),
                      TRUE ~ asked),
    correct = case_when(str_ends(correct, "s") ~ str_sub(correct, 1, -2),
                        TRUE ~ correct)
  ) %>% 
  make_full_confusion_heatmap() +
  labs(#title="Confusion Matrix for Main Occupation",
    caption="Note: Mutiple occupations may have been detected for a response, in which case a single occupations was sampled for probing",
    y="Occupation Detected + Probed", x="Occupation Confirmed") +
  theme(axis.text.x = element_text(angle=60))
ggsave("figures/a_heatmap_occu2.png", height=12, width=12)


## for respondents who were given a confirmation probe, how many 
## said 'no', but then selected that same category in the confirm-Other probe
## (suggests confusion)
conf_other_same <- 
  map_dfr(c("issue","econd_reason","news","occu"),
          function(exp) {
            treat_conf %>%
              rename_at(paste0(exp,"_confirm_other"), ~"confirm_other") %>%
              rename_at(paste0(exp,"_to_confirm"), ~"to_confirm") %>%
              filter(to_confirm != "none_of_the_above") %>%                   
              summarise(pct = mean(confirm_other==to_confirm, na.rm=T),
                        count = sum(confirm_other==to_confirm, na.rm=T)) %>%
              mutate(exp = exp)
          })

# Comparison of Seed Response Quality (Probe vs. No Probe) -----------------
message("Comparison of Seed Response Quality (Probe vs. No Probe)")

if (!exists("coded_qual_pre_post_resolved_with_covariates")) {
  coded_qual_pre_post_resolved_with_covariates <- bind_rows(
    coded_qual_resolved %>%
      mutate(after_probe = treat=="Elaboration", coding_round = 1) %>%
      left_join(combined %>% 
                  select(-matches("(issue|occu|news|econ)"), matches("probe_flag")),
                by = "respondent_id"),
    coded_qual_pre_resolved %>%
      mutate(after_probe = FALSE, coding_round = 2) %>%
      left_join(combined %>% 
                  select(-matches("(issue|occu|news|econ)"), matches("probe_flag")),
                by = "respondent_id")
  )
}

coded_qual_pre_post_resolved_with_covariates %>%
  filter(after_probe == FALSE, treat == "Elaboration") %>%
  mutate(got_probe = case_when(
    q == "econd" & (econd_elab_probe_flag | econd_rel_probe_flag) ~ T,
    q == "issue" & (issue_elab_probe_flag | issue_rel_probe_flag) ~ T,
    q == "news" & (news_elab_probe_flag | news_rel_probe_flag) ~ T,
    q == "occu" & (occu_elab_probe_flag | occu_rel_probe_flag) ~ T,
    
    q == "econd" & (!econd_elab_probe_flag & !econd_rel_probe_flag) ~ F,
    q == "issue" & (!issue_elab_probe_flag & !issue_rel_probe_flag) ~ F,
    q == "news" & (!news_elab_probe_flag & !news_rel_probe_flag) ~ F,
    q == "occu" & (!occu_elab_probe_flag & !occu_rel_probe_flag) ~ F,
    
    q == "econd" & (is.na(econd_elab_probe_flag) & is.na(econd_rel_probe_flag)) ~ F,
    q == "issue" & (is.na(issue_elab_probe_flag) & is.na(issue_rel_probe_flag)) ~ F,
    q == "news" & (is.na(news_elab_probe_flag) & is.na(news_rel_probe_flag)) ~ F,
    q == "occu" & (is.na(occu_elab_probe_flag) & is.na(occu_rel_probe_flag)) ~ F,
    
    .default = NA
  )) %>% #count(q, got_probe)
  group_by(q, got_probe) %>%
  summarise(across(all_of(quality_criteria), ~mean(.x, na.rm=T)), 
            .groups = "drop") %>%
  # focus only on questions where there is a comparison group
  pivot_longer(cols=all_of(quality_criteria)) %>%
  mutate(got_probe = ifelse(got_probe, "Probe Triggered", "No Probe Triggered"),
         name = stringr::str_to_title(name),
         q = exp_ordered_labels_f(q)) %>%
  mutate(name = factor(name, levels=stringr::str_to_title(quality_criteria))) %>%
  ggplot(aes(x=got_probe, y=value)) +
  facet_grid(name ~ q, scales = "free_x", labeller=label_wrap_gen(width = 10)) +
  labs(caption = str_wrap(paste("Note: For occupation experiment,",
                                "probes were administered for every response."), 120),
       x = "Textbot Behavior After Seed Response") +
  geom_bar(stat="identity", fill="grey") +
  geom_text(aes(label=ifelse(
    value != 0, sprintf("%0.1f%%", 100*value), "N/A"
  ),
  x=ifelse(
    value != 0, got_probe, "Probe Triggered"
  ),
  y=value, 
  hjust=case_when(
    q == "Pref. News" & name == "Explanatory" ~ 1.5, 
    q == "Main Occu." & name == "Explanatory" ~ 0.5,
    TRUE ~ 0.5),
  vjust=ifelse(value > 0.5, 1.5, -.25)), 
  color="black", size=4) +
  coord_cartesian(clip = "off") +
  scale_y_continuous(labels = scales::percent_format(1),
                     name = "Rate of Incidence in Seed Response") +
  theme_bw() +
  theme(legend.position = "none",
        panel.background = element_rect(),
        plot.caption = element_text(lineheight = 1, size = 10, hjust = 0),
        strip.text = element_text(size=22, lineheight=1),
        strip.text.y = element_text(angle=0, hjust = 0),
        axis.text.y = element_text(size=14),
        axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
        axis.title = element_text(size=18, lineheight=1),
        panel.spacing = unit(1, "lines"),
        plot.margin = margin(t = 10, r = 10, b = 10, l = 10))

ggsave(here(paste0("figures/a_compare_seed_resp_qual", IMG_FILE_EXT)), 
       height=12, width=12)

# Effects of Elab/Rel Probing on Coded Response Quality (Pre-Post) ---------
message("Effects of Elab/Rel Probing on Coded Response Quality (Pre-Post)")

if (!exists("coded_qual_pre_post_resolved_with_covariates")) {
  coded_qual_pre_post_resolved_with_covariates <- bind_rows(
    coded_qual_resolved %>%
      mutate(after_probe = treat=="Elaboration", coding_round = 1) %>%
      left_join(combined %>% 
                  select(-matches("(issue|occu|news|econ)"), matches("probe_flag")),
                by = "respondent_id"),
    coded_qual_pre_resolved %>%
      mutate(after_probe = FALSE, coding_round = 2) %>%
      left_join(combined %>% 
                  select(-matches("(issue|occu|news|econ)"), matches("probe_flag")),
                by = "respondent_id")
  )
}

if (!exists("coded_qual_all_results")) {
  coded_qual_all_results <- map_dfr(
    quality_criteria,
    .progress = TRUE,
    function(out) {
      reg_subsets <- list("Most Imp.\nIssue"="issue", 
                          "Econ. Cond.\n(Reason)"="econd",
                          "Pooled"=c("issue","econd"))
      lapply(1:length(reg_subsets), function(k) {
        # for some reason, NAs with `lm_robust()` (not enough power?), sticking to `lm()`
        bind_rows(
          # Diff-in-Means: no controls
          coded_qual_pre_post_resolved_with_covariates %>%
            rename_at(out, ~"y") %>%
            filter(q %in% reg_subsets[[k]]) %>%
            mutate(y = as.numeric(y)) %>%
            filter((treat == "Elaboration" & after_probe) | 
                     (treat == "Control" & !after_probe)) %>%
            lm(y ~ treat, data = .) %>%
            tidy() %>%
            mutate(type = "No Controls", design = "Diff-in-Means")
          ,
          # Diff-in-Means: control for covariates
          coded_qual_pre_post_resolved_with_covariates %>%
            rename_at(out, ~"y") %>%
            filter(q %in% reg_subsets[[k]]) %>%
            mutate(y = as.numeric(y)) %>%
            filter((treat == "Elaboration" & after_probe) | 
                     (treat == "Control" & !after_probe)) %>%
            filter(!is.na(after_probe)) %>%
            lm(y ~ treat + work_status + gender + hh_inc5 + educ5 + log(age+1) + device_type, data = .) %>%
            tidy() %>%
            mutate(type = "Controls", design = "Diff-in-Means")
          ,
          # Pre-Post: no controls
          coded_qual_pre_post_resolved_with_covariates %>%
            rename_at(out, ~"y") %>%
            filter(q %in% reg_subsets[[k]]) %>%
            mutate(y = as.numeric(y)) %>%
            filter(!is.na(after_probe)) %>%
            lm(y ~ after_probe, data = .) %>%
            tidy() %>%
            mutate(type = "No Controls", design = "Pre-Post")
          ,
          # Pre-Post: control for covariates
          coded_qual_pre_post_resolved_with_covariates %>%
            rename_at(out, ~"y") %>%
            filter(q %in% reg_subsets[[k]]) %>%
            mutate(y = as.numeric(y)) %>%
            filter(!is.na(after_probe)) %>%
            lm(y ~ after_probe + work_status + gender + hh_inc5 + educ5 + log(age+1) + device_type, data = .) %>%
            tidy() %>%
            mutate(type = "Controls", design = "Pre-Post")
          # ,
          # # Diff-in-Diff: no controls
          # coded_qual_pre_post_resolved_with_covariates %>%
          #   rename_at(out, ~"y") %>%
          #   filter(q %in% reg_subsets[[k]]) %>%
          #   mutate(y = as.numeric(y)) %>%
          #   filter(!is.na(after_probe)) %>%
          #   lm(y ~ after_probe*treat, data = .) %>%
          #   tidy() %>%
          #   mutate(type = "No Controls", design = "Diff-in-Diff")
          # ,
          # # Diff-in-Diff: control for covariates
          # coded_qual_pre_post_resolved_with_covariates %>%
          #   rename_at(out, ~"y") %>%
          #   filter(q %in% reg_subsets[[k]]) %>%
          #   mutate(y = as.numeric(y)) %>%
          #   filter(!is.na(after_probe)) %>%
          #   lm(y ~ after_probe*treat + work_status + gender + hh_inc5 + educ5 + log(age+1) + device_type, data = .) %>%
          #   tidy() %>%
          #   mutate(type = "Controls", design = "Diff-in-Diff")
        ) %>%
          mutate(outcome = out %>% 
                   stringr::str_replace_all("_", " ") %>%
                   stringr::str_to_title(),
                 subset = names(reg_subsets)[[k]])
      })
    })
}

coded_qual_all_results %>%
  filter(design == "Pre-Post") %>%
  mutate(outcome = as_factor(outcome),
         subset = as_factor(subset),
         type = as_factor(type)) %>%
  adj_bhq(alpha = global_stat_sig_level) %>%
  filter(grepl("after_probe", term), !grepl("treat", term)) %>%
  ggplot(aes(x=type, y=estimate, ymin=conf.low, ymax=conf.high,
             group=type, shape=type, alpha=ifelse(stat.sig, stat_sig_label, stat_insig_label))) +
  geom_hline(yintercept = 0, lty = 2, alpha = 0.5) +
  geom_pointrange(size = 1, position = position_dodge(width=0.5)) +
  geom_text(aes(label = sprintf("%.1f%%", estimate*100) %>%
                  ifelse(abs(estimate) < 0.001, "≈0%", .) %>%
                  ifelse(estimate > 0 & abs(estimate) >= 0.001, paste0("+",.), .), 
                y = ifelse(estimate > 0.3, conf.low, conf.high),
                vjust = ifelse(estimate > 0.3, 1.5, -0.5)),
            position = position_dodge(width=0.75), size=5) +
  labs(caption = str_wrap(paste("Note: Higher levels correspond to more positive evaluations.",
                                "Control covariates include work status, gender, education, age, and device type."), 150)) +
  scale_y_continuous(name = "OLS Estimate of Probe in Elab/Rel Condition\n(95% CI with BHq Correction)",
                     labels = scales::percent_format(1)) +
  scale_x_discrete(name = "Estimate", labels = \(.) gsub("Treatment Group","",.)) +
  scale_alpha_manual(values=c(1, 0.2), name = "P-value:") +
  scale_fill_identity(name = "P-value:") +
  scale_shape(name = "Specification:") +
  facet_grid(outcome ~ subset, scales="free_x") +
  theme_bw() +
  theme(strip.text = element_text(size=22, angle=0),
        strip.text.y = element_text(angle=0, hjust=0),
        legend.position = "bottom",
        legend.title = element_text(size=16),
        legend.text = element_text(size=14),
        plot.caption = element_text(lineheight = 1, size = 10, hjust = 0),
        panel.background = element_rect(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_text(size=18),
        axis.title = element_text(size=16, lineheight=1))

ggsave(here(paste0("figures/a_fx_qual_criteria_pp", IMG_FILE_EXT)), 
       height=12, width=12)

# Het FX of Elab/Rel Probing on Human-Coded Quality Criteria --------------
message("Het FX of Elab/Rel Probing on Human-Coded Quality Criteria")

coded_qual_resolved_with_covariates <- coded_qual_resolved %>%
  # reverse code
#  mutate(across(all_of(quality_criteria_rev), ~case_match(.x, 0 ~ 1, 1 ~ 0, .default = NA))) %>%
  left_join(combined %>% 
              select(-matches("(issue|occu|news|econ)"), matches("probe_flag")),
            by = "respondent_id")

coded_qual_hetfx_results <- map_dfr(
  quality_criteria,
  .progress = TRUE,
  function(out) {
    message(out)
    reg_subsets <- list("Most Imp. Issue"="issue", 
                        "Econ. Cond. (Reason)"="econd",
                        "Both"=c("issue","econd"))
    map_dfr(1:length(reg_subsets), 
            function(k) {
              # for some reason, NAs with `lm_robust()` (not enough power?), sticking to `lm()`
              purrr::map_dfr(
                c("work_status","gender","educ5","age5","device_type","hh_inc5"),
                function(var) {
                  # message(" ",var)
                  purrr::map_dfr(
                    unique(sort(combined[[var]])) %>%
                      na.omit() %>%
                      str_filter("Tablet") %>%
                      str_filter("(Other|Refused|IDK)"),
                    function(val) {
                      # message("  ",val)
                      d <- coded_qual_resolved_with_covariates %>%
                        filter_at(var, ~.x == val)%>%
                        rename_at(out, ~"y") %>%
                        filter(q %in% reg_subsets[[k]]) %>%
                        mutate(y = as.numeric(y)) %>%
                        filter(!is.na(y)) %>%
                        filter(treatment %in% c("Control", "Elab/Rel."))

                      d %>%
                        lm(y ~ treatment, data = .) %>%
                        tidy() %>%
                        mutate(type = "No Controls") %>%
                        mutate(outcome = out %>% 
                                 stringr::str_replace_all("_", " ") %>%
                                 stringr::str_to_title(),
                               subset = names(reg_subsets)[[k]]) %>%
                        mutate(val = val,
                               var = var_labels[[var]],
                               # @Soubhik: let's negate the estimates for
                               # het fx regressions when there isn't
                               # any outcome variation (explaining why no/tiny std
                               # errors and a zero estimate) .. see below, lets
                               # put an 'N/A' text as a placeholder 
                               outcome_variation = length(unique(d$y)) > 1,
                               estimate = ifelse(!outcome_variation, NA, estimate),
                               p.value = ifelse(!outcome_variation, NA, p.value),
                               std.error = ifelse(!outcome_variation, NA, std.error))
                    })
                })
            })
  }) 

coded_qual_hetfx_results %>%
  filter(subset == "Most Imp. Issue") %>%
  mutate(outcome = as_factor(outcome %>%
                               gsub("Comprehensible","Comprehen\n-sible", .)),
         subset = as_factor(subset),
         type = as_factor(type)) %>%
  adj_bhq(alpha = global_stat_sig_level) %>%
  filter(grepl("treat", term)) %>%
  # filter(!is.na(p.value)) %>% #@Soubhik: make sure to NOT filter on is.na(p.value)
  ggplot(aes(y=val, x=estimate, xmin=conf.low, xmax=conf.high,
             # group=type, shape=type, 
             color = ifelse(is.na(stat.sig), "NA", ifelse(stat.sig, stat_sig_label, stat_insig_label)),
             fill = ifelse(is.na(stat.sig), "NA", ifelse(stat.sig, stat_sig_label, stat_insig_label)))) +
  geom_hline(yintercept = 0, lty = 2, alpha = 0.5) +
  geom_text(aes(label=ifelse(is.na(estimate), "N/A", ""), x=0, y=val), 
            color="grey") +
  geom_pointrange(size = 1, position = position_dodge(width=0.5)) +
  labs(caption = str_wrap(paste("Note: Higher levels correspond to more positive evaluations.",
                                "Control covariates include work status, gender, education, age, and device type.",
                                # @Soubhik: add the note below to each caption
                                "N/A indicates not a large enough subgroup or enough outcome variation",
                                "to estimate heterogeneity."), 120)) +
  scale_x_continuous(name = "OLS Estimate of Elab/Rel. Treatment Effect\n(95% CI with BHq Correction)",
                     labels = scales::number_format(accuracy=0.01)) +
  geom_vline(xintercept = 0, lty=2, alpha = 0.5) +
  scale_y_discrete(name = "Subgroup", limits = rev) +
  scale_color_manual(
    values = c("p < 0.05" = "black", "p >= 0.05" = "darkgray", "NA" = "grey"), 
    name = "P-value:") +
  scale_fill_manual(
    values = c("p < 0.05" = "black", "p >= 0.05" = "darkgray", "NA" = "grey"), 
    name = "P-value:") +
  guides(
    color = guide_legend(override.aes = list(shape = c(NA, 16, 16)))
    ) +
#  scale_shape(name = "Specification:") +
  facet_grid(var ~ outcome, scales = "free", space = "free_y", labeller=label_wrap_gen(width = 10)) +
  theme_bw() +
  theme(strip.text.x = element_text(size=14, lineheight=1),
        legend.position = "bottom",
        legend.title = element_text(size=12),
        legend.text = element_text(size=12),
        plot.caption = element_text(lineheight = 1, size=10, hjust=0),
        panel.background = element_rect(),
        strip.text.y = element_text(angle=0, size=18, hjust=0.1),
        axis.text.y = element_text(size=12),
        axis.text.x = element_text(size=12, hjust=1, angle = 45),
        axis.title = element_text(size=16, lineheight=1),
        panel.spacing.x = unit(1, "lines"),
        panel.spacing.y = unit(1, "lines"),
        plot.margin = margin(t = 10, r = 10, b = 10, l = 10))

ggsave(here(paste0("figures/a_hetfx_qual_criteria", IMG_FILE_EXT)), 
       height=12, width=12)

# Het FX of Elab/Rel Probing on NLP-Based Informational Measures ----------
message("Het FX of Elab/Rel Probing on NLP-Based Informational Measures")

nlp_measures_l <- nlp_measures %>%
  mutate(response = case_when(
    treatment == "Control" & grepl("econd_reason_1", response) ~ "econd_last_response",
    TRUE ~ response
  )) %>%
  mutate(exp = gsub("_(combined|last)_response", "", response) %>%
           exp_ordered_labels_f()) %>%
  pivot_longer(names_to = "measure",
               values_to = "val",
               unique_words:shannon_entropy) %>% 
  select(-starts_with("f"))

nlp_measures_l_with_covariates <- nlp_measures_l %>%
  left_join(combined %>% 
              select(-matches("(issue|occu|news|econ)"), matches("probe_flag")),
            by = c("treatment","respondent_id"))

table(nlp_measures_l_with_covariates$exp)

nlp_measure_names_main <- c("unique_words"="Unique Words",
                            "total_words"="Total Words",
#                            "lexical_diversity"="Lexical Diversity",
#                            "kl_divergence"="KL Divergence",
#                            "shannon_entropy"="Shannon Entropy",
                            nlp_measure_names[c(3,4,5)])

nlp_hetfx_results <- map_dfr(
  1:length(names(nlp_measure_names_main)),
  .progress = TRUE,
  function(k) {
    message(names(nlp_measure_names_main)[k])
    map_dfr(c("issue", "news", "occu"),
            .progress = TRUE,
            function(e) {
              message(" ", e)
              if (e == "occu")
                vars <- c("gender","educ5","age5","device_type","hh_inc5")
              else
                vars <- c("work_status","gender","educ5","age5","device_type","hh_inc5")
              map_dfr(
                vars,
                function(var) {
                  message(" ",var)
                  map_dfr(
                    unique(sort(combined[[var]])) %>%
                      na.omit() %>%
                      str_filter("Tablet") %>%
                      str_filter("(Other|Refused|IDK)"),
                    function(val) {
                      d <- nlp_measures_l_with_covariates %>%
                        dplyr::filter_at(var, ~.x == val)
                      
                      d_1 <- d %>% 
                        filter(!(grepl("(News|Occu)", exp) & grepl("Control", treatment))) %>%
                        filter(!(grepl("(Issue|Econ)", exp) & grepl("Conf", treatment))) %>%
                        filter(grepl("last_response", response) & !(grepl("why_frustrated", response))) %>%
                        filter(measure == names(nlp_measure_names_main)[k] &
                                 exp == exp_ordered_labels[[e]]) %>%
                        rename(y = val) %>%
                        filter(!is.na(y)) %>%
                        mutate(y = (y - mean(y, na.rm=T)/sd(y, na.rm=T)))
                      d_2 <- d %>% 
                        filter(!(grepl("(News|Occu)", exp) & grepl("Control", treatment))) %>%
                        filter(!(grepl("(Issue|Econ)", exp) & grepl("Conf", treatment))) %>%
                        mutate(response = case_when(
                          treatment == "Control" & grepl("econd_last", response) ~ "econd_combined_response",
                          TRUE ~ response
                        )) %>%
                        filter(grepl("combined_response", response) & !(grepl("why_frustrated", response))) %>%
                        filter(measure == names(nlp_measure_names_main)[k] &
                                 exp == exp_ordered_labels[[e]]) %>%
                        rename(y = val) %>%
                        filter(!is.na(y)) %>%
                        mutate(y = (y - mean(y, na.rm=T)/sd(y, na.rm=T)))
                      bind_rows(
                        if (all(is.finite(d_1$y))) {
                          d_1 %>%
                            lm(y ~ treatment, data = .) %>%
                            tidy(conf.level = 0.95, conf.int = TRUE) %>%
                            mutate(type = "Post Probe",
                                   outcome_variation = length(unique(d_1$y)) > 1)
                        } else { data.frame() },
                        if (all(is.finite(d_2$y))) {
                          d_2 %>%
                            lm(y ~ treatment, data = .) %>%
                            tidy(conf.level = 0.95, conf.int = TRUE) %>%
                            mutate(type = "Combined",
                                   outcome_variation = length(unique(d_2$y)) > 1)
                        } else { data.frame() }
                      ) %>%
                        mutate(measure = nlp_measure_names_main[[k]] %>% 
                                 stringr::str_wrap(width = 10),
                               val = val,
                               var = var_labels[[var]],
                               exp = e,
                               estimate = ifelse(!outcome_variation, NA, estimate),
                               p.value = ifelse(!outcome_variation, NA, p.value),
                               std.error = ifelse(!outcome_variation, NA, std.error))
                    })
                })
            })
  })


nlp_hetfx_results %>%
  filter(exp == "issue", type == "Post Probe") %>%
  mutate(exp = exp_ordered_labels_f(exp),
         val = as_factor(val),
         measure = as_factor(measure)) %>%
  mutate(stat.sig = p.value < global_stat_sig_level) %>%
  filter(grepl("treatment", term)) %>%
  adj_bhq(alpha = global_stat_sig_level) %>%
  filter(grepl("treat", term)) %>%
  filter(!(grepl("other", val, ignore.case=TRUE))) %>%
  filter(!(grepl("tablet", val, ignore.case=TRUE))) %>%
  ggplot(aes(y=val, x=estimate, xmin=conf.low, xmax=conf.high,
             # shape=type,
             # group=type,
#             color=case_when(
#               conf.low > 0 ~ "green3",
#               !stat.sig ~ "black",
#               conf.high < 0 ~ "red"
#             ),
             # group=type, shape=type,
#             alpha=ifelse(stat.sig, stat_sig_label, stat_insig_label))) 
              color = ifelse(is.na(stat.sig), "NA", ifelse(stat.sig, stat_sig_label, stat_insig_label)),
              fill = ifelse(is.na(stat.sig), "NA", ifelse(stat.sig, stat_sig_label, stat_insig_label))))+
  scale_x_continuous(name = "OLS Estimate of Elab/Rel Group Membership on NLP Measure\n(95% CI with BHq Correction)",
                     labels = scales::number_format(accuracy = 0.01)) +
  scale_y_discrete(name = "Subgroup", limits = rev) +
  labs(caption = str_wrap(paste("Note: Outcomes are standardized via z-transformation.",
                                "N/A indicates not a large enough subgroup or enough outcome variation",
                                "to estimate heterogeneity."), 120)) +
  geom_vline(xintercept = 0, lty = 2, alpha = 0.5) +
  geom_pointrange(size = 1, position = position_dodge(0.5)) +
  geom_text(aes(label=ifelse(is.na(estimate), "N/A", ""), x=0, y=val),
            color = "grey") +
  facet_grid(var ~ measure, scales = "free", space = "free_y", labeller=label_wrap_gen(width = 10)) +
  scale_color_manual(values = c("p < 0.05" = "black", "p >= 0.05" = "darkgray", "NA" = "grey"), name = "P-value:") +
  scale_fill_manual(values = c("p < 0.05" = "black", "p >= 0.05" = "darkgray", "NA" = "grey"), name = "P-value:") +
  # scale_shape(name = "Measure:") +
  # scale_alpha_manual(values=c(0.2), name = "P-value:") +
  # scale_fill_identity(name = "Effect:") +
  # scale_color_identity(name = "Effect:") +
  # scale_shape(name = "Outcome\nMeasure:") +
  theme_bw() + 
  theme(strip.text = element_text(size=18, lineheight=1),
        legend.position = "bottom",
        legend.title = element_text(size=12),
        legend.text = element_text(size=12),
        plot.caption = element_text(lineheight = 1, size=10, hjust=0),
        panel.background = element_rect(),
        strip.text.y = element_text(angle=0, hjust=0.1),
        axis.text.y = element_text(size=12),
        axis.text.x = element_text(size=12, hjust=1, angle = 45),
        axis.title = element_text(size=16, lineheight=1),
        panel.spacing.x = unit(1, "lines"),
        panel.spacing.y = unit(1, "lines"),
        plot.margin = margin(t = 10, r = 10, b = 10, l = 10))

ggsave(here(paste0("figures/a_hetfx_nlp_measures", IMG_FILE_EXT)), 
       height=12, width=12)

# Het FX of Probing on Respondent Attrition -------------------------------
# message("Het FX of Probing on Respondent Attrition")
# 
# # Note: not running because no real heterogeneity to speak of
#
# ## No real heterogeneity here
# nr_types <- list("Dropout (Right After)"="_dropout_rightafter",
#                  "Dropout (After)"="_dropout_after")
# 
# map_dfr(
#   c("issue","econd_reason","news","occu"),
#   .progress = TRUE,
#   function(exp) {
#     message(exp)
#     map_dfr(
#       c("work_status","gender","educ5","age5","device_type","hh_inc5"),
#       function(var) {
#         # message(" ",var)
#         map_dfr(
#           unique(sort(combined[[var]])) %>%
#             na.omit() %>%
#             str_filter("Tablet") %>%
#             str_filter("(Other|Refused|IDK)"),
#           function(val) {
#             # message("  ",val)
#             # just diff-in-means
#             combined %>%
#               filter_at(var, ~.x == val) %>%
#               rename_at(paste0(exp, "_dropout_after"), ~"y") %>%
#               filter(treatment %in% c("Control", "Elab/Rel.")) %>%
#               filter(!is.na(y)) %>%
#               lm(y ~ treatment, data = .) %>%
#               tidy() %>%
#               # View() %>% 
#               mutate(val = val,
#                      var = var_labels[[var]],
#                      exp = exp,
#                      outcome_variation = length(unique(d$y)) > 1)
#           })
#       }) %>%
#       adj_bhq(alpha = global_stat_sig_level)
#   }) %>%
#   mutate_at(c("outcome","var"), ~.x %>%
#               stringr::str_replace("_"," ") %>%
#               stringr::str_replace("age5","age") %>%
#               stringr::str_replace("work status","employed") %>%
#               stringr::str_replace("type","")) %>%
#   mutate(val = as_factor(val),
#          exp = exp_ordered_labels_f(exp)) %>%
#   filter(grepl("treat", term)) %>%
#   filter(!(grepl("other", val, ignore.case=TRUE))) %>%
#   filter(!(grepl("tablet", val, ignore.case=TRUE))) %>%
# #  View() %>% 
#   ggplot(aes(y=val, x=estimate, xmin = conf.low, xmax = conf.high,
#              # group=type, shape=type, 
#              color = ifelse(is.na(stat.sig), "NA", ifelse(stat.sig, stat_sig_label, stat_insig_label)),
#              fill = ifelse(is.na(stat.sig), "NA", ifelse(stat.sig, stat_sig_label, stat_insig_label)))) +
#   geom_pointrange(size = 0.5, position = position_dodge(width=0.5), shape = 21) +
#   # geom_text(aes(label = sprintf("%0.2f", estimate), 
#   #               x = ifelse(estimate > 0, conf.low, conf.high),
#   #               hjust = ifelse(estimate > 0, 1.5, -1.25)), 
#   #           position = position_dodge(width=0.75), size=5) +
#   labs(caption = str_wrap(paste("Note: Higher levels correspond to more positive evaluations.",
#                                 "Control covariates include work status, gender, education, age, and device type.",
#                                 "Outcomes are all normalized to the [0-1] range.",
#                                 "Excluded subgroups in 'Other' categories and who answered by tablet due to small sample sizes."), 120)) +
#   scale_x_continuous(name = "Difference Between Elab/Rel. and Control Group Conditions\n(OLS Coefficient with Corrected 95% CI)",
#                      labels = scales::number_format(accuracy = 0.01)) +
#   scale_y_discrete(name = "Subgroup", limits = rev) +
#   #scale_alpha_manual(values=c(1, 1, 0.2), name = "P-value:") +
# #  scale_fill_identity(name = "P-value:") +
#   scale_color_manual(values = c("p < 0.05" = "black", "p >= 0.05" = "darkgray", "NA" = "darkgray"), name = "P-value:") +
#   scale_fill_manual(values = c("p < 0.05" = "black", "p >= 0.05" = "darkgray", "NA" = "white"), name = "P-value:") +
#   guides(color = guide_legend(override.aes = list(fill = c("white", "black", "darkgray")))) +
#   geom_vline(xintercept = 0, lty=2, alpha = 0.5) +
#   facet_grid(var ~ exp, scales = "free", space = "free_y", labeller=label_wrap_gen(width = 10)) +
#   theme_bw() +
#   theme(strip.text = element_text(size=22, lineheight=1),
#         legend.position = "bottom",
#         legend.title = element_text(size=12),
#         legend.text = element_text(size=12),
#         plot.caption = element_text(lineheight = 1, size=10, hjust=0),
#         strip.text.y = element_text(angle=0, hjust = 0.1),
#         panel.background = element_rect(),
#         axis.text.y = element_text(size=16),
#         axis.text.x = element_text(size=12, angle=45, hjust=1),
#         axis.title = element_text(size=18, lineheight=1),
#         panel.spacing.x = unit(1, "lines"),
#         panel.spacing.y = unit(1, "lines"),
#         plot.margin = margin(t = 10, r = 10, b = 10, l = 10))
# 
# # Note: not saving because no real heterogeneity to speak of
# # ggsave(here(paste0("figures/a_hetfx_attrition", IMG_FILE_EXT)), 
# #        height=12, width=12)


# Het FX of Probing on Self-Reported User Experience ----------------------
message("Het FX of Probing on Self-Reported User Experience")

map_dfr(c("quality","ease","frustration","satisfaction"),
        .progress = T,
        function(out) {
          message(out)
          map_dfr(c("work_status","gender","educ5","age5","device_type","hh_inc5"),
                  function(var) {
                    # message(" ",var)
                    map_dfr(unique(sort(combined[[var]])) %>%
                              na.omit() %>%
                              str_filter("Tablet") %>%
                              str_filter("(Other|Refused|IDK)"),
                            function(val) {
                              # message("  ",val)
                              # just diff-in-means
                              combined %>%
                                filter(treatment %in% c("Control", "Elab/Rel.")) %>%
                                rename_at(out, ~"y") %>%
                                mutate(y = as.numeric(y)) %>%
                                filter(!is.na(y)) %>%
                                mutate(y = y/max(y, na.rm=T)) %>%
                                filter_at(var, ~.x == val) %>%
                                lm(y ~ treatment, data = .) %>%
                                tidy() %>%
                                mutate(val = val,
                                       var = var_labels[[var]],
                                       outcome = out)
                            })
                  }) %>%
            adj_bhq(alpha = global_stat_sig_level)
        }) %>%
  mutate_at(c("outcome","var"), ~.x %>%
              stringr::str_replace("_"," ") %>%
              stringr::str_replace("age5","age") %>%
              stringr::str_replace("work status","employed") %>%
              stringr::str_replace("type","")) %>%
  mutate(val = as_factor(val),
         outcome = stringr::str_to_title(outcome)) %>%
  filter(grepl("treat", term)) %>%
  filter(!(grepl("other", val, ignore.case=TRUE))) %>%
  filter(!(grepl("tablet", val, ignore.case=TRUE))) %>%
  ggplot(aes(y=val, x=estimate, xmin=conf.low, xmax=conf.high,
             # group=type, shape=type, 
             color=ifelse(stat.sig, stat_sig_label, stat_insig_label))) +
  geom_pointrange(size = 1, position = position_dodge(width=0.5)) +
  # geom_text(aes(label = sprintf("%0.2f", estimate), 
  #               x = ifelse(estimate > 0, conf.low, conf.high),
  #               hjust = ifelse(estimate > 0, 1.5, -1.25)), 
  #           position = position_dodge(width=0.75), size=5) +
  labs(caption = str_wrap(paste("Note: Higher levels correspond to more positive evaluations.",
                                "Control covariates include work status, gender, education, age, and device type.",
                                "Outcomes are all normalized to the [0-1] range.",
                                "Excluded subgroups in 'Other' categories and who answered by tablet due to small sample sizes."), 120)) +
  scale_x_continuous(name = "Difference Between Elab/Rel. and Control Group Conditions\n(OLS Coefficient with Corrected 95% CI)",
                     labels = scales::number_format(accuracy = 0.01)) +
  scale_y_discrete(name = "Subgroup", limits = rev) +
  # scale_alpha_manual(values=c(1, 0.2, 0), name = "P-value:") +
  scale_color_manual(values = c("p < 0.05" = "black", "p >= 0.05" = "darkgray", "NA" = "darkgray"), name = "P-value:") +
  scale_fill_manual(values = c("p < 0.05" = "black", "p >= 0.05" = "darkgray", "NA" = "white"), name = "P-value:") +
  geom_vline(xintercept = 0, lty=2, alpha = 0.5) +
  facet_grid(var ~ outcome, scales = "free", space = "free_y", labeller=label_wrap_gen(width=10)) +
  theme_bw() +
  theme(strip.text = element_text(size=22, lineheight=1),
        legend.position = "bottom",
        legend.title = element_text(size=12),
        legend.text = element_text(size=12),
        plot.caption = element_text(lineheight = 1, size = 10, hjust=0),
        strip.text.y = element_text(angle=0, hjust = 0),
        panel.background = element_rect(),
        axis.text.y = element_text(size=16),
        axis.text.x = element_text(size=12, angle=45, hjust=1),
        axis.title = element_text(size=18, lineheight=1),
        panel.spacing.x = unit(1, "lines"),
        panel.spacing.y = unit(1, "lines"),
        plot.margin = margin(t = 10, r = 10, b = 10, l = 10))

ggsave(here(paste0("figures/a_hetfx_user_experience", IMG_FILE_EXT)), 
       height=12, width=12)


message("DONE.")
